ndarrayのfortran flagと実際のメモリ配列について¶

要約¶

ndarrayで.Tで転置を行った場合、実際のメモリ上の配列は変更されず、.flagsで確認できるC形式配列(row-major)であるかFortran形式配列(column-major)であるかのフラグが入れ替わるのみ
np.saveで保存する場合も上記フラグは保たれたまま保存されるため、実際のレイアウトを転置させたい場合はndarrayに対し.copy()を読んだ上で保存する必要がある

実験¶

import numpy as np

2x2のndarrayを生成する。

a = np.array([[1,2], [3,4]])
a

array([[1, 2],
       [3, 4]])

.Tで転置を取ると行と列が入れ替わる。

a.T

array([[1, 3],
       [2, 4]])

元のndarrayの.flagsを表示すると、C_CONTIGUOUSがTrueとなっている一方、F_CONTIGUOUSがFalseとなっていることが分かる。C_CONTIGUOUSはC形式配列であるかを表すフラグであり、F_CONTIGUOUSはFortran形式配列であるかを表すフラグである。

a.flags

  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

一方.Tを呼んだ後のndarrayの.flagsを表示すると、C_CONTIGUOUSがFalseに変化し、一方F_CONTIGUOUSがTrueとなっている。

a.T.flags

  C_CONTIGUOUS : False
  F_CONTIGUOUS : True
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

`np.save`での保存時の形式¶

.Tで転置を取ると.flagsのC形式配列であるかFortran形式配列であるかのフラグが入れ替わることが分かったが、実際の格納形式がどうなっているかを知るために、np.saveを用いて保存した.npyファイルに対して!hexdump -Cvを実行することでバイナリの中身を調べる。ここでhexdumpコマンドの-CはASCII表示を行うためのオプションであり、-vは全データを表示するためのオプションである(おそらくBSD系のhexdumpのみ)。

まず、np.saveで転置前後のndarrayを保存する。

np.save("a.npy", a)
np.save("a.T.npy", a.T)

これらをnp.loadで再度読み込み直すと、.flagsは読み込んだ後も保持されていることが分かる。

np.load("a.npy").flags

  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

np.load("a.T.npy").flags

  C_CONTIGUOUS : False
  F_CONTIGUOUS : True
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

ディスクへダンプした際の実際の格納形式がどうなっているかをhexdumpコマンドを用いて調べる。

!hexdump -Cv a.npy

00000000  93 4e 55 4d 50 59 01 00  76 00 7b 27 64 65 73 63  |.NUMPY..v.{'desc|
00000010  72 27 3a 20 27 3c 69 38  27 2c 20 27 66 6f 72 74  |r': '<i8', 'fort|
00000020  72 61 6e 5f 6f 72 64 65  72 27 3a 20 46 61 6c 73  |ran_order': Fals|
00000030  65 2c 20 27 73 68 61 70  65 27 3a 20 28 32 2c 20  |e, 'shape': (2, |
00000040  32 29 2c 20 7d 20 20 20  20 20 20 20 20 20 20 20  |2), }           |
00000050  20 20 20 20 20 20 20 20  20 20 20 20 20 20 20 20  |                |
00000060  20 20 20 20 20 20 20 20  20 20 20 20 20 20 20 20  |                |
00000070  20 20 20 20 20 20 20 20  20 20 20 20 20 20 20 0a  |               .|
00000080  01 00 00 00 00 00 00 00  02 00 00 00 00 00 00 00  |................|
00000090  03 00 00 00 00 00 00 00  04 00 00 00 00 00 00 00  |................|
000000a0

!hexdump -Cv a.T.npy

00000000  93 4e 55 4d 50 59 01 00  76 00 7b 27 64 65 73 63  |.NUMPY..v.{'desc|
00000010  72 27 3a 20 27 3c 69 38  27 2c 20 27 66 6f 72 74  |r': '<i8', 'fort|
00000020  72 61 6e 5f 6f 72 64 65  72 27 3a 20 54 72 75 65  |ran_order': True|
00000030  2c 20 27 73 68 61 70 65  27 3a 20 28 32 2c 20 32  |, 'shape': (2, 2|
00000040  29 2c 20 7d 20 20 20 20  20 20 20 20 20 20 20 20  |), }            |
00000050  20 20 20 20 20 20 20 20  20 20 20 20 20 20 20 20  |                |
00000060  20 20 20 20 20 20 20 20  20 20 20 20 20 20 20 20  |                |
00000070  20 20 20 20 20 20 20 20  20 20 20 20 20 20 20 0a  |               .|
00000080  01 00 00 00 00 00 00 00  02 00 00 00 00 00 00 00  |................|
00000090  03 00 00 00 00 00 00 00  04 00 00 00 00 00 00 00  |................|
000000a0

ヘッダーのfortran_orderが変化していることが分かる。更にbashのプロセス置換で両者のdiffを見てみる。

%%bash
diff <(hexdump -Cv a.npy) <(hexdump -Cv a.T.npy)

3,5c3,5
< 00000020  72 61 6e 5f 6f 72 64 65  72 27 3a 20 46 61 6c 73  |ran_order': Fals|
< 00000030  65 2c 20 27 73 68 61 70  65 27 3a 20 28 32 2c 20  |e, 'shape': (2, |
< 00000040  32 29 2c 20 7d 20 20 20  20 20 20 20 20 20 20 20  |2), }           |
---
> 00000020  72 61 6e 5f 6f 72 64 65  72 27 3a 20 54 72 75 65  |ran_order': True|
> 00000030  2c 20 27 73 68 61 70 65  27 3a 20 28 32 2c 20 32  |, 'shape': (2, 2|
> 00000040  29 2c 20 7d 20 20 20 20  20 20 20 20 20 20 20 20  |), }            |

ヘッダーのみ差分が存在し、実際のデータが格納されている後半部は変化していないことが分かる。

転置した状態のまま保存する方法¶

ndarrayに対して.copy()を呼んだ上で保存すれば良い。下記に実例を示す。まず、.copy()を呼んだ後の.flagsを表示する。

a.T.copy().flags

  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

転置したにもかかわらずC_CONTIGUOUSがTrueになっていることが分かる。更にこのndarrayをnp.saveで保存し、hexdumpで中身を見てみる。

np.save("a.T.copy.npy", a.T.copy())

!hexdump -C a.T.copy.npy

00000000  93 4e 55 4d 50 59 01 00  76 00 7b 27 64 65 73 63  |.NUMPY..v.{'desc|
00000010  72 27 3a 20 27 3c 69 38  27 2c 20 27 66 6f 72 74  |r': '<i8', 'fort|
00000020  72 61 6e 5f 6f 72 64 65  72 27 3a 20 46 61 6c 73  |ran_order': Fals|
00000030  65 2c 20 27 73 68 61 70  65 27 3a 20 28 32 2c 20  |e, 'shape': (2, |
00000040  32 29 2c 20 7d 20 20 20  20 20 20 20 20 20 20 20  |2), }           |
00000050  20 20 20 20 20 20 20 20  20 20 20 20 20 20 20 20  |                |
*
00000070  20 20 20 20 20 20 20 20  20 20 20 20 20 20 20 0a  |               .|
00000080  01 00 00 00 00 00 00 00  03 00 00 00 00 00 00 00  |................|
00000090  02 00 00 00 00 00 00 00  04 00 00 00 00 00 00 00  |................|
000000a0

C_CONTIGUOUSがTrue、F_CONTIGUOUSがFalseになっていたことからも分かるように、ヘッダーのfortran_orderがFalseになっている。

当然、再度np.loadで読み込んだ場合でもこの.flagsは保たれている。

np.load("a.T.copy.npy").flags

  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

元のndarrayと、転置した上で.copy()を呼んだ後のndarrayをnp.saveしたバイナリを比較してみる。

%%bash
diff <(hexdump -Cv a.npy) <(hexdump -Cv a.T.copy.npy)

9,10c9,10
< 00000080  01 00 00 00 00 00 00 00  02 00 00 00 00 00 00 00  |................|
< 00000090  03 00 00 00 00 00 00 00  04 00 00 00 00 00 00 00  |................|
---
> 00000080  01 00 00 00 00 00 00 00  03 00 00 00 00 00 00 00  |................|
> 00000090  02 00 00 00 00 00 00 00  04 00 00 00 00 00 00 00  |................|

両者ともにC形式配列となっているため、純粋にデータ部の並びのみが異なる。

同様に転置を取った場合と、その上で.copy()を呼んだ場合のバイナリを比較してみる。

%%bash
diff <(hexdump -Cv a.T.npy) <(hexdump -Cv a.T.copy.npy)

3,5c3,5
< 00000020  72 61 6e 5f 6f 72 64 65  72 27 3a 20 54 72 75 65  |ran_order': True|
< 00000030  2c 20 27 73 68 61 70 65  27 3a 20 28 32 2c 20 32  |, 'shape': (2, 2|
< 00000040  29 2c 20 7d 20 20 20 20  20 20 20 20 20 20 20 20  |), }            |
---
> 00000020  72 61 6e 5f 6f 72 64 65  72 27 3a 20 46 61 6c 73  |ran_order': Fals|
> 00000030  65 2c 20 27 73 68 61 70  65 27 3a 20 28 32 2c 20  |e, 'shape': (2, |
> 00000040  32 29 2c 20 7d 20 20 20  20 20 20 20 20 20 20 20  |2), }           |
9,10c9,10
< 00000080  01 00 00 00 00 00 00 00  02 00 00 00 00 00 00 00  |................|
< 00000090  03 00 00 00 00 00 00 00  04 00 00 00 00 00 00 00  |................|
---
> 00000080  01 00 00 00 00 00 00 00  03 00 00 00 00 00 00 00  |................|
> 00000090  02 00 00 00 00 00 00 00  04 00 00 00 00 00 00 00  |................|

この場合、ヘッダーもデータ部も異なっていることが分かる。

なぜこの記事を書いたか¶

とある.npyファイルのデータ部分を直接読み込むCプログラムを触っていた際に、.Tで転置した.npyを読み込んでも読み込んだデータ部分が全く変化していないということが発生したためである。当該プログラムではヘッダーのfortran_orderを完全に無視していたため、.Tしただけのndarrayを保存した.npyではデータ部分が転置されておらず、意図しない動作となってしまっていた。これを解決するためには前述のように.copy()を呼んだ上で保存した.npyファイルを読み込む必要があった。

2018/5/19にPyCon mini Osakaで「librosaで始める音楽情報検索」というタイトルで発表しました¶

公開が遅くなりましたが先月開催されたPyCon mini Osakaで発表した内容を公開します。この記事もAzure Notebooksで既に公開済ではありましたが、nikolaのipynb公開機能を使ってブログポストとして投稿しておきたいと思います。

librosaで始める音楽情報検索¶

大橋宏正(twitter: @wrist)

この資料¶

Azure Notebooksで作成
- https://bit.ly/2rRjYSj

自己紹介¶

大橋宏正(@wrist)
- 某メーカで音響信号処理の業務に従事
- 2013年ごろから大阪PRML読書会を主催(現在休止中)
- PyData Osakaオーガナイザ

PyCon mini Osakaのテーマ¶

「Pythonでなにしよう？」
- Webアプリ
- 科学技術計算
Pythonで音楽の分析を行いたい！

この発表¶

Pythonで音楽情報検索(MIR)を行うためのライブラリlibrosaの紹介

音楽情報検索(MIR)とは？¶

Music Information Retrievalの略
音楽情報に関する処理全般
- ビートトラッキング、楽器音分離、音楽推薦など
コミュニティ
- ismir.net
- 音楽情報科学研究会(www.sigmus.jp)

MIRをPythonで行う¶

Pythonで音楽データの処理が必要
読み書き
- wave, audioread, pysoundfile, scipy.io, ...
線形代数、信号処理、フーリエ変換
- numpy, scipy(signal, fftpack, ...)
可視化
- matplotlib, seaborn
これらを組み合わせることでさまざまな処理を実現

参考¶

音響信号処理についてはQiitaに過去に諸々記事を書いてます

librosa¶

音楽情報処理を行うためのPythonライブラリ
- コロンビア大学にかつて存在したLabROSAで開発
- 2/17に0.6.0がリリース
github
- https://github.com/librosa/librosa
ドキュメント
- https://librosa.github.io/librosa/

librosaを学ぶための資料¶

Scipy 2015 conf.の資料、動画
- http://conference.scipy.org/proceedings/scipy2015/brian_mcfee.html
ギャラリー
- https://github.com/librosa/librosa_gallery
- https://librosa.github.io/librosa_gallery/
チュートリアル
- https://github.com/librosa/tutorial
- https://librosa.github.io/librosa/tutorial.html

librosaの特徴¶

MIRに必要となる典型的な処理が多く実装
低い学習障壁
一貫性、後方互換性、モジュール性

# pipでインストール
!pip install librosa

Requirement already satisfied: librosa in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (0.6.0)
Requirement already satisfied: scipy>=0.14.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (1.1.0)
Requirement already satisfied: resampy>=0.2.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (0.2.0)
Requirement already satisfied: numpy>=1.8.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (1.14.2)
Requirement already satisfied: decorator>=3.0.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (4.3.0)
Requirement already satisfied: audioread>=2.0.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (2.1.5)
Requirement already satisfied: six>=1.3 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (1.11.0)
Requirement already satisfied: scikit-learn!=0.19.0,>=0.14.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (0.19.1)
Requirement already satisfied: joblib>=0.7.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (0.11)
Requirement already satisfied: numba>=0.32 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from resampy>=0.2.0->librosa) (0.38.0)
Requirement already satisfied: llvmlite>=0.23.0dev0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from numba>=0.32->resampy>=0.2.0->librosa) (0.23.0)

librosaの具体的な紹介¶

適宜デモを交えながらlibrosaを紹介
発表の流れ
- 音データの短時間周波数分析
- 特徴量抽出
- オンセット検出、ビートトラッキング

# 本発表では以下を仮定
import numpy as np
import scipy as sp
import scipy.signal as sg
import scipy.fftpack as fft
import matplotlib as mpl
import matplotlib.pyplot as plt
import IPython.display

# version
print("numpy version", np.version.full_version)
print("scipy version", sp.version.full_version)
print("matplotlib version", mpl.__version__)

numpy version 1.14.2
scipy version 1.1.0
matplotlib version 2.2.2

# librosa以下の各モジュールをimport
import librosa
import librosa.display
import librosa.feature
import librosa.effects
import librosa.onset
import librosa.beat

サウンドデータの読み込み¶

librosa.load
- ファイルパスを与えることでデータを読み込み
- デフォルトでモノラル化、サンプリングレート22050Hzとなる
- バックエンドにaudioreadライブラリを使用
librosa.display.waveplotで波形データの可視化

# デモ音源のロード＆再生
# y: モノラル音声データ, sr: サンプリングレート
y, sr = librosa.load(librosa.util.example_audio_file())
IPython.display.Audio(data=y, rate=sr)

# 波形の描画
# 縦軸レンジを自動的に指定し、横軸を秒数表示にしてくれる
librosa.display.waveplot(y, sr)

<matplotlib.collections.PolyCollection at 0x115239668>

サウンドデータのデジタル表現¶

現実の音はアナログデータ
- 時間方向にも振幅方向にも連続
コンピュータで音を取り扱うためには離散化が必要
- 時間方向: サンプリング
- 振幅方向: 量子化

b_pt = int(10.0*sr); e_pt = int((10.0+0.002)*sr)
fig = plt.figure(1, figsize=(14,6)); ax = fig.add_subplot(111)
ax.plot(y[b_pt:e_pt]); ax.plot(y[b_pt:e_pt], "o"); ax.grid();
ax.set_xlabel("time [pt]"); ax.set_ylabel("amplitude")

Text(0,0.5,'amplitude')

# サンプリングレート: 1秒辺りのサンプル数
print("sampling rate: {0}".format(sr))
print("length: {0}[pt]={1}[s]".format(y.shape[0], y.shape[0]/sr))

sampling rate: 22050
length: 1355168[pt]=61.45886621315193[s]

音声データの周波数分析¶

音声データの長さはデータにより異なる
固定長の短時間波形に切り出した上で処理を行う
短時間波形に対し離散フーリエ変換を適用することで周波数分析を行う

# strideトリックを使った信号を短時間波形フレームの配列に変換する関数
from numpy.lib.stride_tricks import as_strided

def split_frames(xs, frame_width=1024, noverlap=2):

    nsample = np.alen(xs)
    nshift = (frame_width // noverlap)
    nframe = int(nsample / nshift) - noverlap + 1
    
    nbyte = xs.dtype.itemsize

    # as_strided(xs, ys_shape, new_stride)
    xs_split = as_strided(xs, (nframe, frame_width), (nbyte * nshift, nbyte))
    
    return xs_split

def plot_overlap_waves():
    import matplotlib.gridspec as gridspec

    frame_width = 2048
    half_width = int(frame_width // 2)
    nframes = 20

    fig = plt.figure(1, figsize=(12,8))
    gs = gridspec.GridSpec(3, nframes + 1)
    
    ax_wav = fig.add_subplot(gs[0, :])
    ax_wav.plot(np.arange(half_width*(nframes+1))/sr, 
                y[:half_width*(nframes+1)])
    ymax = np.max(y[:half_width*(nframes+1)])
    ax_wav.set_xlim([0, half_width*(nframes+1)/sr])
    ax_wav.set_ylim([-ymax, ymax]); ax_wav.grid()

    win = sg.get_window("hanning", int(frame_width))
    xs = np.arange(frame_width)

    for i, y_cut in enumerate(split_frames(y, int(frame_width), 2)[:nframes]):
        y_win = win * y_cut
        ax_wav.plot(np.arange(half_width*i, half_width*i + frame_width)/sr, 0.2 * win)
        ax = None
        if i % 2 == 0:
            ax = fig.add_subplot(gs[1, i:i+2])
            ax.plot(xs[:half_width], y_win[:half_width])
            ax.plot(xs[half_width:], y_win[half_width:])
            ax.set_xlim([0, frame_width])
        else:
            ax = fig.add_subplot(gs[2, i:i+2])
            ax.plot(xs[half_width:], y_win[half_width:])
            ax.plot(xs[:half_width], y_win[:half_width])
            ax.set_xlim([0, frame_width])

        ax.tick_params(labelleft="off",left="off") # y軸の削除
        ax.set_ylim([-ymax, ymax]); ax.grid()

    fig.subplots_adjust(wspace=0.0)

# 周波数分析時には通常窓関数を乗算した短時間波形を使用
# 切り出しはここでは1/2区間をオーバーラップさせている
plot_overlap_waves()

/Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: Passing one of 'on', 'true', 'off', 'false' as a boolean is deprecated; use an actual boolean (True/False) instead.
  warnings.warn(message, mplDeprecation, stacklevel=1)

用語の整理¶

sr = sampling rate [Hz = pt/s]
- 22050 [Hz]
frame = 各短時間波形
n_fft = 各フレームに含まれるサンプル数
- 2048 [pt]
hop_length = フレーム間のサンプル数
- 512 [pt]

def plot_overlap_freqs():
    import matplotlib.gridspec as gridspec

    frame_width = 2048
    half_width = int(frame_width // 2)
    nframes = 20

    fig = plt.figure(1, figsize=(12,8))
    gs = gridspec.GridSpec(3, nframes + 1)
    
    ax_wav = fig.add_subplot(gs[0, :])
    ax_wav.plot(y[:half_width*(nframes+1)])
    ymax = np.max(y[:half_width*(nframes+1)])
    ax_wav.set_xlim([0, half_width*(nframes+1)])
    ax_wav.set_ylim([-ymax, ymax]); ax_wav.grid()

    win = sg.get_window("hanning", int(frame_width))
    xs = np.arange(frame_width)
    w = np.linspace(0.0, 2.0 * np.pi, frame_width, endpoint=False)
    freqs = w[:frame_width//2+1] * (sr / 2) / np.pi
    for i, y_cut in enumerate(split_frames(y, int(frame_width), 2)[:20]):
        y_win = win * y_cut
        Y_abs = np.abs(fft.fft(y_win, frame_width)[:frame_width//2+1])
        ax_wav.plot(np.arange(half_width*i, half_width*i + frame_width), 0.2 * win)
        ax = None
        if i % 2 == 0:
            ax = fig.add_subplot(gs[1, i:i+2])
            ax.semilogy(Y_abs, freqs)
        else:
            ax = fig.add_subplot(gs[2, i:i+2])
            ax.semilogy(Y_abs, freqs)

        if i > 1:
            ax.tick_params(labelleft="off",left="off") # y軸の削除
        else:
            ax.tick_params(labelbottom="off",bottom="off") # x軸の削除
            ax.set_ylabel("frequency [Hz]")

        ax.set_xlabel("amplitude")
        ax.set_ylim([10, sr/2])
        ax.grid()

    fig.subplots_adjust(wspace=0.0)

# 各短時間波形を周波数分析したグラフを描画
# 縦軸が周波数、横軸が振幅スペクトル
plot_overlap_freqs()

/Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: Passing one of 'on', 'true', 'off', 'false' as a boolean is deprecated; use an actual boolean (True/False) instead.
  warnings.warn(message, mplDeprecation, stacklevel=1)

def plot_spectrogram():
    frame_width = 2048
    half_width = int(frame_width // 2)
    nframes = 20

    fig = plt.figure(1, figsize=(12,8))
    
    ax_wav = fig.add_subplot(211)
    ax_wav.plot(np.arange(half_width*(nframes+1))/sr,
                y[:half_width*(nframes+1)])
    ymax = np.max(y[:half_width*(nframes+1)])
    ax_wav.set_xlim([0, half_width*(nframes+1)/sr])
    ax_wav.set_ylim([-ymax, ymax]); ax_wav.grid()

    win = sg.get_window("hanning", int(frame_width))
    xs = np.arange(frame_width)
    w = np.linspace(0.0, 2.0 * np.pi, frame_width, endpoint=False)
    freqs = w[:frame_width//2+1] * (sr / 2) / np.pi

    Zs = []
    for i, y_cut in enumerate(split_frames(y, int(frame_width), 2)[:20]):
        y_win = win * y_cut
        Y_abs = np.abs(fft.fft(y_win, frame_width)[:frame_width//2+1])
        ax_wav.plot(np.arange(half_width*i, half_width*i + frame_width), 0.2 * win)
        Zs.append(20.0*np.log10(Y_abs))
    
    ax = fig.add_subplot(212)
    Xs = np.arange((nframes+1)) / sr / 2
    pcol = ax.pcolormesh(Xs, freqs, np.array(Zs).T, cmap='magma')
    cax = fig.colorbar(pcol)
    cax.set_clim(-80, 20)

    fig.subplots_adjust(wspace=0.0)

# 短時間波形を周波数分析したものをスペクトログラムとして表示
plot_spectrogram()

librosaにおける短時間周波数分析¶

短時間離散フーリエ変換(Short time fourie transform)
- librosa.stft
Constant-Q Transform(CQT)
- librosa.cqt
分析結果の可視化
- librosa.specshow

# librosaでは短時間波形分析をメソッド一発で実行可能
D = librosa.stft(y)

# D[k, l]
# k: 周波数ビン  l: フレームインデックス
# kは実信号に対する2048ptフーリエ変換の結果なので1025ptのみ保持
D.shape

(1025, 2647)

離散フーリエ変換$D[k,l]$の算出法¶

$$D[k, l]=\sum_{n=0}^{{\textrm N_{\rm fft}-1}} w[n]\cdot x\left[\left(l\cdot {\rm M}-\frac{N_{\rm fft}}{2}\right)+n\right]\exp\left\{-j\frac{2\pi k n}{N_{\rm fft}}\right\}$$

* $x[n]$は入力信号、$w[n]$は窓関数
* 周波数ビン$k$、フレームインデックス$l$、$M$はhop length
* $lM$を中心に前後$N/2$点に渡り離散フーリエ変換を実行

fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

log_power = librosa.amplitude_to_db(np.abs(D), ref=np.max)
librosa.display.specshow(log_power, x_axis="time", y_axis="linear")
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x115a2aba8>

# 縦軸を対数表示
fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

log_power = librosa.amplitude_to_db(np.abs(D), ref=np.max)
librosa.display.specshow(log_power, x_axis="time", y_axis="log")
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x1c3344cef0>

Constant-Q変換¶

$$C[k]= \frac{1}{N[k]} \sum_{n=0}^{{\textrm N[k]-1}} w[n, k]\cdot x\left[n\right]\exp\left\{-j\frac{2\pi Q[k] n}{N[k]}\right\}$$

特徴
- 分析する周波数間隔を対数的に配置
- 低域は細かく分析し、高域は荒く分析
- 音楽信号の分析に適する

C = librosa.cqt(y, sr)

# STFTと比較して縦軸の周波数データ数は圧倒的に少ない
# STFTは周波数間隔が等幅なのに対してCQTだと対数間隔となるため
print(D.shape)
print(C.shape)

(1025, 2647)
(84, 2647)

fmin = 32.70  # C1に相当するHz
fmax = fmin * ((2**(1/12)) ** 84)
print(fmax)
print(librosa.hz_to_note(fmin))
print(librosa.hz_to_note(fmax))

4185.600000000015
C1
C8

fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

log_cqt_power = librosa.amplitude_to_db(np.abs(C), ref=np.max)
librosa.display.specshow(log_cqt_power, x_axis="time", y_axis="cqt_hz")
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x115933cf8>

fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

log_cqt_power = librosa.amplitude_to_db(np.abs(C), ref=np.max)
librosa.display.specshow(log_cqt_power, x_axis="time", y_axis="cqt_note")
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x11589d518>

fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

# デフォルトだと80dBまでだがtop_dbを付けるとレンジが縮まる
log_cqt_power = librosa.amplitude_to_db(np.abs(C), ref=np.max, top_db=40)
librosa.display.specshow(log_cqt_power, x_axis="time", y_axis="cqt_note")
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x1c1f823518>

クロマベクトル¶

12音階それぞれの強さを要素に持つベクトル
CQTで算出した対数周波数帯域ごとの強度を元に計算
librosa.feature.chroma_cqt

chroma = librosa.feature.chroma_cqt(C=C, sr=sr)

fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

librosa.display.specshow(chroma, x_axis="time", y_axis="chroma")
plt.colorbar()

/Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages/librosa/display.py:657: UserWarning: Trying to display complex-valued input. Showing magnitude instead.
  warnings.warn('Trying to display complex-valued input. '

<matplotlib.colorbar.Colorbar at 0x1159a3128>

# STFTの結果Dからクロマベクトルを算出
chroma_stft = librosa.feature.chroma_stft(S=D, sr=sr)

fig = plt.figure(1, figsize=(12,4))
ax = fig.add_subplot(2,1,1)
librosa.display.specshow(chroma_stft, y_axis='chroma')
plt.title('chroma_stft'); plt.colorbar()
ax= fig.add_subplot(2,1,2)
librosa.display.specshow(chroma, y_axis='chroma', x_axis='time')
plt.title('chroma_cqt'); plt.colorbar()
plt.tight_layout()

/Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages/librosa/display.py:657: UserWarning: Trying to display complex-valued input. Showing magnitude instead.
  warnings.warn('Trying to display complex-valued input. '

その他特徴量抽出¶

メルスペクトル、MFCC、tonnetzなども算出可能
$\Delta$特徴量の計算も可能

M = librosa.feature.melspectrogram(y=y, sr=sr)
mfcc = librosa.feature.mfcc(y=y, sr=sr)
tonnetz = librosa.feature.tonnetz(y=y, sr=sr)

`librosa.effects`モジュール¶

音源分離
- librosa.hpss
  - 調波成分と打音成分に分解
時間伸縮、ピッチシフト
- librosa.time_stretch
- librosa.pitch_shift

y_harmonic, y_percussive = librosa.effects.hpss(y)

IPython.display.Audio(y, rate=sr)

IPython.display.Audio(y_harmonic, rate=sr)

IPython.display.Audio(y_percussive, rate=sr)

# それぞれに対するCQTを求める
C      = librosa.cqt(y=y, sr=sr)
C_harm = librosa.cqt(y=y_harmonic, sr=sr)
C_perc = librosa.cqt(y=y_percussive, sr=sr)

fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

log_cqt_power= librosa.amplitude_to_db(np.abs(C), ref=np.max)
log_cqt_harm = librosa.amplitude_to_db(np.abs(C_harm), ref=np.max(np.abs(C)))
log_cqt_perc = librosa.amplitude_to_db(np.abs(C_perc), ref=np.max(np.abs(C)))

ax = fig.add_subplot(3,1,1); librosa.display.specshow(log_cqt_power, x_axis="time", y_axis="cqt_hz")
ax = fig.add_subplot(3,1,2); librosa.display.specshow(log_cqt_harm , x_axis="time", y_axis="cqt_hz")
ax = fig.add_subplot(3,1,3); librosa.display.specshow(log_cqt_perc , x_axis="time", y_axis="cqt_hz")

<matplotlib.axes._subplots.AxesSubplot at 0x1139ae2e8>

# 2倍の速度に時間伸縮
y_fast = librosa.effects.time_stretch(y, 2.0)
IPython.display.Audio(y_fast, rate=sr)

# 3度上にピッチシフト
y_third = librosa.effects.pitch_shift(y, sr, n_steps=4)
IPython.display.Audio(y_third, rate=sr)

# 減五度にピッチシフト
y_tritone = librosa.effects.pitch_shift(y, sr, n_steps=-6)
IPython.display.Audio(y_tritone, rate=sr)

onset検出とビートトラッキング¶

librosa.onsetモジュール、librosa.beatモジュールを使用
- 出音の強度をonset.onset_strengthで算出
- 強度を元にonset箇所をonset.onset_detectで検出
- 強度を元にbeat.beat_trackでテンポとビートを推定

onset_envelope = librosa.onset.onset_strength(y, sr)
onsets = librosa.onset.onset_detect(onset_envelope=onset_envelope)

fig = plt.figure(1, figsize=(12,8))
ax = fig.add_subplot(2, 1, 1)
ax.plot(onset_envelope, label="onset strength")
ax.vlines(onsets, 0, np.max(onset_envelope), color="red", alpha=0.25, label="onsets")
ax.set_xlim([0, len(onset_envelope)]); ax.legend(loc=1)
ax = fig.add_subplot(2, 1, 2)
librosa.display.waveplot(y, sr)

<matplotlib.collections.PolyCollection at 0x113aa6c88>

tempo, beats = librosa.beat.beat_track(onset_envelope=onset_envelope)

print(tempo)

129.19921875

fig = plt.figure(1, figsize=(12,8))
ax = fig.add_subplot(1, 1, 1)
ax.plot(onset_envelope, label="onset strength")
ax.vlines(onsets, 0, 0.8 * np.max(onset_envelope), color="red", alpha=0.5, linestyle="dotted", label="onsets")
ax.vlines(beats, 0, np.max(onset_envelope), color="black", alpha=0.25, label="beats")
ax.set_xlim([0, len(onset_envelope)]); ax.legend(loc=1)

<matplotlib.legend.Legend at 0x114fc1128>

クリック音を生成し同時に再生¶

MIRの評価に用いられるmir_evalライブラリを使用

!pip install mir_eval

Requirement already satisfied: mir_eval in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (0.4)
Requirement already satisfied: six in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from mir_eval) (1.11.0)
Requirement already satisfied: numpy>=1.7.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from mir_eval) (1.14.2)
Requirement already satisfied: scipy>=0.9.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from mir_eval) (1.1.0)
Requirement already satisfied: future in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from mir_eval) (0.16.0)

# http://colinraffel.com/publications/ismir2014mir_eval.pdf
import mir_eval

beat_times = librosa.frames_to_time(beats)
y_click = mir_eval.sonify.clicks(beat_times, sr, length=len(y))

IPython.display.Audio(y_click, rate=sr)

IPython.display.Audio(y + y_click, rate=sr)

特徴量のビート同期¶

librosa.util.syncを用いるとビートとその他特徴量を同期させることができる
- 各ビートの間の成分を平均化したり中央値を取ったりできる

tempo, beat_frames = librosa.beat.beat_track(y=y_percussive, sr=sr)
# 13次元のMFCCおよびΔを抽出
mfcc = librosa.feature.mfcc(y=y, sr=sr, n_mfcc=13)
mfcc_delta = librosa.feature.delta(mfcc)

print(beat_frames.shape)
print(mfcc.shape)
print(mfcc_delta.shape)

(98,)
(13, 2647)
(13, 2647)

# ビート間の平均値を用いてMFCCとΔMFCCをビートに同期
beat_mfcc_delta = librosa.util.sync(np.vstack([mfcc, mfcc_delta]),
                                    beat_frames)

beat_mfcc_delta.shape

(26, 99)

# 調波成分からクロマベクトルを推定
chromagram = librosa.feature.chroma_cqt(y=y_harmonic, sr=sr)

# ビート間の中央値を用いてクロマベクトルをビートに同期
beat_chroma = librosa.util.sync(chromagram, beat_frames, aggregate=np.median)

# ビート同期特徴量(クロマベクトル、MFCC、ΔMFCC)を結合
beat_features = np.vstack([beat_chroma, beat_mfcc_delta])

# MFCC 13dim, ΔMFCC 13dim, クロマベクトル(12音階; 12dim)
beat_features.shape

(38, 99)

fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(211)
librosa.display.specshow(log_cqt_power, x_axis="time", y_axis="cqt_hz")
log_beat_features = librosa.amplitude_to_db(np.abs(beat_features), ref=np.max)
ax = fig.add_subplot(212)
ax.pcolormesh(log_beat_features)

<matplotlib.collections.QuadMesh at 0x113a0ac50>

まとめ¶

librosaを用いた音楽情報検索(MIR)について説明
- 特徴量抽出、可視化、音源分離、ビート推定
音楽データに対する探索的分析の手段の一つ

以下補足およびおまけ¶

`librosa.decompose`による音源分離¶

NMFによるモノラル音源分離を行う
- スペクトログラムを行列と見てscikit-learnのNMFに投げる
- アクティベーション行列と基底行列に分解
- 教師なし音源分離が可能

NMF¶

非負値行列$\textbf{X}$を$\textbf{X} \simeq \textbf{WH} $と分解する手法
- Hを基底行列、Wをアクティベーションと呼ぶ
scikit-learnのNMFでは下記目的関数を最小化

$$0.5 * ||\textbf{X} - \textbf{WH}||_{Fro}^2

\alpha l1_{ratio} ||vec(W)||_1
\alpha l1_{ratio} ||vec(H)||_1
0.5 \alpha (1 - l1{ratio}) ||W||{Fro}^2
0.5 \alpha (1 - l1{ratio}) ||H||{Fro}^2 $$
通常最適化には補助関数法が使われる

D = librosa.stft(y)
S, phase = librosa.magphase(D)
comps, acts = librosa.decompose.decompose(S, n_components=8, sort=True)
S_approx = comps @ acts

print(comps.shape)
print(acts.shape)
print(S.shape, S_approx.shape)

(1025, 8)
(8, 2647)
(1025, 2647) (1025, 2647)

def plot_decomposed(S, comps, acts, S_approx):
    fig = plt.figure(figsize=(10,8))
    ax = fig.add_subplot(311)
    librosa.display.specshow(librosa.amplitude_to_db(S, ref=np.max), y_axis='log', x_axis='time')
    ax.set_title('Input spectrogram'); plt.colorbar(format='%+2.0f dB')

    ax = fig.add_subplot(323)
    librosa.display.specshow(librosa.amplitude_to_db(comps, ref=np.max), y_axis='log')
    plt.colorbar(format='%+2.0f dB'); ax.set_title('Components')

    ax = fig.add_subplot(324)
    librosa.display.specshow(acts, x_axis='time')
    ax.set_ylabel('Components'); ax.set_title('Activations')
    plt.colorbar()

    ax = fig.add_subplot(313)
    librosa.display.specshow(librosa.amplitude_to_db(S_approx, ref=np.max), y_axis='log', x_axis='time')
    plt.colorbar(format='%+2.0f dB')
    ax.set_title('Reconstructed spectrogram')

    fig.tight_layout()

plot_decomposed(S, comps, acts, S_approx)

S_partly = comps[:, :3] @ acts[:3, :]

plot_decomposed(S, comps[:, :3], acts[:3, :], S_partly)

y_partly = librosa.istft(S_partly * phase)
IPython.display.Audio(y_partly, rate=sr)

decomposeの用途¶

視聴用途の音源分離には(このままでは)向かない
- 振幅スペクトルのみを分解するが、位相情報は元信号の元を使うため歪が大きい
- 基底数を増加した上で適切なクラスタリングを施して使用したり、事前学習した楽器別の基底などを用いて分離を行うなどの工夫が必要
音響イベント検出のための特徴量に用いることも可能